#include <sfv/array1d.h>
#include <sfv/csr_matrix.h>
#include <sfv/detail/format_utils.h>

#include <stack>
#include <deque>

namespace sfv
{
	template< typename MatrixType, typename Array1, typename Array2 >
	int node_betweenness(const MatrixType& graph, Array1& res, const Array2& vids, int cutoff ) 
	{
		typedef typename MatrixType::index_type IndexType;
		typedef typename MatrixType::value_type ValueType;

		long int no_of_nodes = graph.num_rows;
		long int no_of_edges = graph.num_entries;

		std::deque<long int> q;

		sfv::array1d<IndexType> row_indices(graph.num_entries);
		sfv::detail::offsets_to_indices(graph.row_offsets, row_indices);

		/* here we go */

		for (IndexType source=0; source<no_of_nodes; source++) 
		{
			sfv::array1d<long int> distance( no_of_nodes, 0 );
			sfv::array1d<unsigned long long int> nrgeo( no_of_nodes, 0 );
			sfv::array1d<double> tmpscore( no_of_nodes, 0 );
			std::stack<long int> stack;

			q.push_back(source);
			nrgeo[source]=1;
			distance[source]=0;

			while (!q.empty()) 
			{
				long int actnode = q.front();
				q.pop_front();

				if (cutoff > 0 && distance[actnode] >= cutoff) continue;

				size_t neino = graph.row_offsets[actnode+1] - graph.row_offsets[actnode];
				sfv::array1d<IndexType> neis(neino);
				IndexType base = graph.row_offsets[actnode];
				for( size_t i = 0; i < neino; i++ )
					neis[i] = graph.column_indices[base+i];

				for (size_t j=0; j < neino; j++) {
					long int neighbor = neis[j];
					if (nrgeo[neighbor] != 0) {
						/* we've already seen this node, another shortest path? */
						if (distance[neighbor]==distance[actnode]+1) {
							nrgeo[neighbor]+=nrgeo[actnode];
						}
					} else {
						/* we haven't seen this node yet */
						nrgeo[neighbor]+=nrgeo[actnode];
						distance[neighbor]=distance[actnode]+1;
						q.push_back(neighbor);
						stack.push(neighbor);
					}
				}
			} /* while !igraph_dqueue_empty */

			/* Ok, we've the distance of each node and also the number of
			shortest paths to them. Now we do an inverse search, starting
			with the farthest nodes. */
			while (!stack.empty()) 
			{
				long int actnode = stack.top();
				stack.pop();
				if (distance[actnode]<=1) { continue; } /* skip source node */

				/* set the temporary score of the friends */
				size_t neino = graph.row_offsets[actnode+1] - graph.row_offsets[actnode];
				sfv::array1d<IndexType> neis(neino);
				IndexType base = graph.row_offsets[actnode];
				for( size_t i = 0; i < neino; i++ )
					neis[i] = graph.column_indices[base+i];

				for (size_t j=0; j<neis.size(); j++) {
					long int neighbor = neis[j];
					if (distance[neighbor]==distance[actnode]-1 && nrgeo[neighbor] != 0) {
						tmpscore[neighbor] += 
							(tmpscore[actnode]+1)*nrgeo[neighbor]/nrgeo[actnode];
					}
				}
			}

			/* Ok, we've the scores for this source */
			for (size_t k=0; k < vids.size(); k++) 
			{
				long int node = vids[k];
				res[k] += tmpscore[node];
				tmpscore[node] = 0.0; /* in case a node is in vids multiple times */
			}

		} /* for source < no_of_nodes */

		/* divide by 2 for undirected graph */
		for (size_t j=0; j< res.size(); j++) 
		{
			res[j] /= ValueType(2.0);
		}

		return 0;
	}

	template< typename MatrixType, typename Array >
	int edge_betweenness(const MatrixType& graph, Array& result, int cutoff) 
	{
		typedef typename MatrixType::index_type IndexType;
		typedef typename MatrixType::value_type ValueType;

		long int no_of_nodes = graph.num_rows;
		long int no_of_edges = graph.num_entries;

		std::deque<long int> q;

		result.resize(no_of_edges);

		sfv::array1d<IndexType> row_indices(graph.num_entries);
		sfv::detail::offsets_to_indices(graph.row_offsets, row_indices);

		/* here we go */

		for (IndexType source=0; source<no_of_nodes; source++) 
		{
			sfv::array1d<long int> distance( no_of_nodes, 0 );
			sfv::array1d<unsigned long long int> nrgeo( no_of_nodes, 0 );
			sfv::array1d<double> tmpscore( no_of_nodes, 0 );
			std::stack<long int> stack;

			q.push_back(source);

			nrgeo[source]=1;
			distance[source]=0;

			while (!q.empty()) 
			{

				long int actnode = q.front();
				q.pop_front();

				if (cutoff > 0 && distance[actnode] >= cutoff ) continue;

				size_t neino = graph.row_offsets[actnode+1] - graph.row_offsets[actnode];
				sfv::array1d<IndexType> neip(neino);
				for( size_t i = 0; i < neino; i++ )
					neip[i] = graph.row_offsets[actnode] + i;

				for (size_t i=0; i<neino; i++) 
				{
					int edge = neip[i], from, to;
					long int neighbor;

					from = row_indices[edge];
					to   = graph.column_indices[edge];
					neighbor = (actnode != from) ? from : to;

					if (nrgeo[neighbor] != 0) {
						/* we've already seen this node, another shortest path? */
						if (distance[neighbor]==distance[actnode]+1) {
							nrgeo[neighbor]+=nrgeo[actnode];
						}
					} else {
						/* we haven't seen this node yet */
						nrgeo[neighbor]+=nrgeo[actnode];
						distance[neighbor]=distance[actnode]+1;
						q.push_back(neighbor);
						stack.push(neighbor);
					}
				}
			} /* while !igraph_dqueue_empty */

			/* Ok, we've the distance of each node and also the number of
			shortest paths to them. Now we do an inverse search, starting
			with the farthest nodes. */
			while (!stack.empty()) 
			{
				long int actnode = stack.top();
				stack.pop();
				if (distance[actnode]<1) { continue; } /* skip source node */

				/* set the temporary score of the friends */
				size_t neino = graph.row_offsets[actnode+1] - graph.row_offsets[actnode];
				sfv::array1d<IndexType> neip(neino);
				for( size_t i = 0; i < neino; i++ )
					neip[i] = graph.row_offsets[actnode] + i;

				for (size_t i=0; i<neino; i++) 
				{
					int from, to;
					long int neighbor;
					long int edgeno = neip[i];

					from = row_indices[edgeno];
					to   = graph.column_indices[edgeno];
					neighbor = (actnode != from) ? from : to;

					if (distance[neighbor]==distance[actnode]-1 &&
						nrgeo[neighbor] != 0) 
					{
						tmpscore[neighbor] +=
							(tmpscore[actnode]+1)*nrgeo[neighbor]/nrgeo[actnode];
						result[edgeno] +=
							(tmpscore[actnode]+1)*nrgeo[neighbor]/nrgeo[actnode];
					}
				}
			}
			/* Ok, we've the scores for this source */
		} /* for source <= no_of_nodes */

		/* divide by 2 for undirected graph */
		for (size_t j=0; j<result.size(); j++) {
			result[j] /= ValueType(2.0);
		}

		return 0;
	}

} // end namespace sfv
